#===========================================================================================
# Article:  Heterogeneity in General Multinomial Choice Models
# Journal:  Statistical Methods & Applications 32: 129–148. doi:10.1007/s10260-022-00642-5 
# Author:   Ingrid Mauerer, Gerhard Tutz
#===========================================================================================

#-----------------------------------------------------------------------------------
# This file contains the commands to run the models and create the tables 
# R version 4.2.2 
# Platform: x86_64-w64-mingw32/x64 (64-bit)
#-----------------------------------------------------------------------------------

#==================
# Packages
#==================
library("mlogit")
library("matlib")
library("dplyr")
library("xtable")
library("car")
library("lmtest")
library("gmnl")

#==================
# Load Data
#==================
load("Data.RDS")

#==============
# Run Program 
#==============
# (for details, see Section A in the Supplementary Material)
source("GHMNL.R")

#=============================================================================================================
# GHMNL(dat, namesglob, Indglob, namescats, Indcats, nameshet, nameresp, k, penalty)
#
# - dat:        data set in long format (k rows for each observation) 
#               the response is indicated by a 1 in the corresponding row, 0 in all other rows
#
# - namesglob:  names of global variables (chooser-specific attributes)
# - Indglob >0: global variables are included, otherwise ignored
# 
# - namescats:  names of category-specific  variables (choice-specific attributes),
#               if no category-specific variables are in the data set, choose any variable and set Indcats=0
# - Indcats >0: category-specific variables are included, otherwise ignored
#
# - nameshet:   names of variables in the heterogeneity term (can be empty)
#               if nameshet <- c()  no heterogeneity term in the model
#
# - namesresp:  name of response variable (vector of 0-1 dummy variables, based on nominal response)
# - k:          number of categories in the response
#
# - pen:        penalty on estimates, choose small, pen = 0.0000001 or even pen=0 (unpenalized) 
#               number >= 0, if penalty=0 pure maximum likelihood
#               if penalty>0 ridge penalized maximum likelihood 
#==============================================================================================================



#==================================================================
# GHMNL Model 1: Platform Divergence
#==================================================================
  
k <- 6   
Indglob <- 0
Indcats <- 1
pen <- 0.000001


#+++++++++++++++++++++++++
# Immigration issue (imm)
#+++++++++++++++++++++++++ 
namesglob <- c("interest","know","edubin") 
namescats <- c("imm")
nameshet  <- c("pimm")
nameresp  <- c("chosen")

imm1 <- GHMNL(dat_imm, namesglob,Indglob,namescats,Indcats,nameshet,nameresp,k,pen)
imm1

#-----------------------
# Get Parameters: 
#-----------------------

#---------------
# global var. 
#---------------
# t-value:
imm1globt<-imm1$parglob/imm1$parglobstd

Imm1_glob <- cbind(imm1$parglob,imm1$parglobstd,imm1globt) 
colnames(Imm1_glob)<-c("coef.", "s.e.", "t-val.")
rownames(Imm1_glob)<-c("SPD", "FDP", "Greens", "Left", "AfD")
Imm1_glob


#---------------
# cats var.
#---------------
# t-value:
imm1catst<-imm1$parcats/imm1$parcatsstd

Imm1_cats<- cbind(imm1$parcats,imm1$parcatsstd,imm1catst) 
colnames(Imm1_cats)<-c("coef.", "s.e.", "t-val.")
rownames(Imm1_cats)<-c("Issue Proximity")
Imm1_cats


#---------------
# het var.
#---------------
# t-value:
imm1hett<-imm1$parhet/imm1$parhetstd

Imm1_het_t<- rbind(imm1$parhet,imm1$parhetstd,imm1hett) 
Imm1_het<-t(Imm1_het_t)
colnames(Imm1_het)<-c("coef.", "s.e.", "t-val.")
rownames(Imm1_het)<-c("Platform Divergence")
Imm1_het


#---------------
# Loglikelihood
#---------------
Imm1_logl<-matrix(NA, 1,3)
rownames(Imm1_logl)<-c("Log-likelihood")
Imm1_logl["Log-likelihood",1]<- -imm1$Loglik
Imm1_logl


#---------------------------
# put the matrices together
#---------------------------
IssueImm<-matrix(NA, 1,3)
rownames(IssueImm)<-c("Immigration Issue")

Constant<-matrix(NA,1,3)
rownames(Constant)<-c("Constants")

Heterogeneity<-matrix(NA,1,3)
rownames(Heterogeneity)<-c("Heterogeneity Term")

Location<-matrix(NA,1,3)
rownames(Location)<-c("Location Term")

outputImm1 <-round(rbind(IssueImm, Location, Imm1_glob, Imm1_cats,Heterogeneity, Imm1_het, Imm1_logl), digits=3)
outputImm1


#+++++++++++++++++++
# Tax issue (tax)
#+++++++++++++++++++
namesglob <- c("interest","know","edubin") 
namescats <- c("tax")
nameshet  <- c("ptax")
nameresp  <- c("chosen")

tax1 <- GHMNL(dat_tax, namesglob,Indglob,namescats,Indcats,nameshet,nameresp,k,pen)
tax1

#-----------------------
# Get Parameters: 
#-----------------------

#---------------
# global var. 
#---------------
# t-value:
tax1globt<-tax1$parglob/tax1$parglobstd

Tax1_glob <- cbind(tax1$parglob,tax1$parglobstd,tax1globt) 
colnames(Tax1_glob)<-c("coef.", "s.e.", "t-val.")
rownames(Tax1_glob)<-c("SPD", "FDP", "Greens", "Left", "AfD")
Tax1_glob


#---------------
# cats var.
#---------------
# t-value:
tax1catst<-tax1$parcats/tax1$parcatsstd

Tax1_cats<- cbind(tax1$parcats,tax1$parcatsstd,tax1catst) 
colnames(Tax1_cats)<-c("coef.", "s.e.", "t-val.")
rownames(Tax1_cats)<-c("Issue Proximity")
Tax1_cats


#---------------
# het var.
#---------------
# t-value:
tax1hett<-tax1$parhet/tax1$parhetstd

Tax1_het_t<- rbind(tax1$parhet,tax1$parhetstd,tax1hett) 
Tax1_het<-t(Tax1_het_t)
colnames(Tax1_het)<-c("coef.", "s.e.", "t-val.")
rownames(Tax1_het)<-c("Platform Divergence")
Tax1_het

#---------------
# Loglikelihood
#---------------
Tax1_logl<-matrix(NA, 1,3)
rownames(Tax1_logl)<-c("Log-likelihood")
Tax1_logl["Log-likelihood",1]<- -tax1$Loglik
Tax1_logl


#---------------------------
# put the matrices together
#---------------------------
IssueTax<-matrix(NA, 1,3)
rownames(IssueTax)<-c("Tax Issue")

outputTax1 <-round(rbind(IssueTax,Location, Tax1_glob, Tax1_cats,Heterogeneity, Tax1_het,Tax1_logl), digits=3)
outputTax1

#++++++++++++++++++++++++++++
# Climate Change issue (cc)
#++++++++++++++++++++++++++++
namesglob <- c("interest","know","edubin") 
namescats <- c("cc")
nameshet  <- c("pcc")
nameresp  <- c("chosen")

cc1 <- GHMNL(dat_cc, namesglob,Indglob,namescats,Indcats,nameshet,nameresp,k,pen)
cc1

#-----------------------
# Get Parameters
#-----------------------

#---------------
# global var. 
#---------------
# t-value:
cc1globt<-cc1$parglob/cc1$parglobstd

Cc1_glob <- cbind(cc1$parglob,cc1$parglobstd,cc1globt) 
colnames(Cc1_glob)<-c("coef.", "s.e.", "t-val.")
rownames(Cc1_glob)<-c("SPD", "FDP", "Greens", "Left", "AfD")
Cc1_glob


#---------------
# cats var.
#---------------
# t-value:
cc1catst<-cc1$parcats/cc1$parcatsstd

Cc1_cats<- cbind(cc1$parcats,cc1$parcatsstd,cc1catst) 
colnames(Cc1_cats)<-c("coef.", "s.e.", "t-val.")
rownames(Cc1_cats)<-c("Issue Proximity")
Cc1_cats


#---------------
# het var.
#---------------
# t-value:
cc1hett<-cc1$parhet/cc1$parhetstd

Cc1_het_t<- rbind(cc1$parhet,cc1$parhetstd,cc1hett) 
Cc1_het<-t(Cc1_het_t)
colnames(Cc1_het)<-c("coef.", "s.e.", "t-val.")
rownames(Cc1_het)<-c("Platform Divergence")
Cc1_het

#---------------
# Loglikelihood
#---------------
Cc1_logl<-matrix(NA, 1,3)
rownames(Cc1_logl)<-c("Log-likelihood")
Cc1_logl["Log-likelihood",1]<- -cc1$Loglik
Cc1_logl


#---------------------------
# put the matrices together
#---------------------------
IssueCC<-matrix(NA, 1,3)
rownames(IssueCC)<-c("Climate Change Issue")

outputCc1 <-round(rbind(IssueCC, Location, Cc1_glob, Cc1_cats,Heterogeneity, Cc1_het,Cc1_logl), digits=3)
outputCc1


#-----------------------------------------------------
# Table 2: Platform Divergence: GHMNL Model Estimates
#-----------------------------------------------------
Table2 <-rbind(outputImm1, outputTax1, outputCc1)
Table2

#==================================================================
# GHMNL Model 2: Full Voter Choice Model
#==================================================================
k <- 6   
Indglob <- 1
Indcats <- 1
pen <- 0.000001

namesglob <- c("kath", "gender", "union")
namescats <- c("tax")
nameshet  <- c("edubin", "gender", "west")
nameresp  <- c("chosen")

tax2 <- GHMNL(dat_tax, namesglob,Indglob,namescats,Indcats,nameshet,nameresp,k,pen)
tax2


#-----------------------
# Get Parameters: 
#-----------------------

#---------------
# global var. 
#---------------
# t-value:
tax2globt<-tax2$parglob/tax2$parglobstd


# Constants
tax2_glob1 <- cbind(tax2$parglob[,1],tax2$parglobstd[,1],tax2globt[,1]) 
colnames(tax2_glob1)<-c("coef.", "s.e.", "t-val.")
rownames(tax2_glob1)<-c("SPD", "FDP", "Greens", "Left", "AfD")
tax2_glob1

# kath 
tax2_glob2 <- cbind(tax2$parglob[,2],tax2$parglobstd[,2],tax2globt[,2]) 
colnames(tax2_glob2)<-c("coef.", "s.e.", "t-val.")
rownames(tax2_glob2)<-c("SPD", "FDP", "Greens", "Left", "AfD")
tax2_glob2

# gender 
tax2_glob3 <- cbind(tax2$parglob[,3],tax2$parglobstd[,3],tax2globt[,3]) 
colnames(tax2_glob3)<-c("coef.", "s.e.", "t-val.")
rownames(tax2_glob3)<-c("SPD", "FDP", "Greens", "Left", "AfD")
tax2_glob3

# union membership 
tax2_glob4 <- cbind(tax2$parglob[,4],tax2$parglobstd[,4],tax2globt[,4]) 
colnames(tax2_glob4)<-c("coef.", "s.e.", "t-val.")
rownames(tax2_glob4)<-c("SPD", "FDP", "Greens", "Left", "AfD")
tax2_glob4


#---------------
# cats var.
#---------------
# t-value:
tax2catst<-tax2$parcats/tax2$parcatsstd

Tax2_cats<- cbind(tax2$parcats,tax2$parcatsstd,tax2catst) 
colnames(Tax2_cats)<-c("coef.", "s.e.", "t-val.")
rownames(Tax2_cats)<-c("Issue Proximity")
Tax2_cats


#---------------
# het var.
#---------------
# t-value:
tax2hett<-tax2$parhet/tax2$parhetstd

Tax2_het_t<- rbind(tax2$parhet,tax2$parhetstd,tax2hett) 
Tax2_het<-t(Tax2_het_t)
colnames(Tax2_het)<-c("coef.", "s.e.", "t-val.")
rownames(Tax2_het)<-c("Education", "Gender (1 female, 0 male)", "Region (1 former West Germany, 0 former East Germany)")
Tax2_het

#---------------
# Loglikelihood
#---------------
Tax2_logl<-matrix(NA, 1,3)
rownames(Tax2_logl)<-c("Log-likelihood")
Tax2_logl["Log-likelihood",1]<- -tax2$Loglik
Tax2_logl

#---------------------------
# put the matrices together
#---------------------------
Cath<-matrix(NA,1,3)
rownames(Cath)<-c("Religious Denomination (1 Catholic, 0 otherwise)")

Gender<-matrix(NA,1,3)
rownames(Gender)<-c("Gender (1 female, 0 male)")

Union<-matrix(NA, 1,3)
rownames(Union)<-c("Union Membership (1 union member, 0 otherwise)")

#----------------------------------------------------------
# Table 3: Full Voter Choice Model: GHMNL Model Estimates
#----------------------------------------------------------

Table3 <-round(rbind(Location, Tax2_cats,Constant, tax2_glob1, Cath, tax2_glob2, 
                   Gender, tax2_glob3, Union, tax2_glob4, 
                   Heterogeneity, Tax2_het,Tax2_logl), digits=3)
Table3


#==================================================================
# Standard MNL Models
#==================================================================
# Model 1:

#+++++++++++++++++++++++++
# Immigration issue (imm)
#+++++++++++++++++++++++++  
immMNL <-mlogit(chosen ~ imm  |1| 0, data=dat_imm, shape="long", alt.var="party", reflevel = "1")
summary(immMNL)

#+++++++++++++++++++
# Tax issue (tax)
#+++++++++++++++++++
taxMNL <-mlogit(chosen ~ tax  |1| 0, data=dat_tax, shape="long", alt.var="party", reflevel = "1")
summary(taxMNL)

#++++++++++++++++++++++++++++
# Climate Change issue (cc)
#++++++++++++++++++++++++++++
ccMNL <-mlogit(chosen ~ cc  |1| 0, data=dat_cc, shape="long", alt.var="party", reflevel = "1")
summary(ccMNL)


# Model 2: 
taxMNL2<-mlogit(chosen ~  tax  | kath + gender + union| 0, data=dat_tax, shape="long", alt.var="party", reflevel = "1")
summary(taxMNL2)



#################################################
# LR tests: Compare GHMNL and Standard MNL model
#################################################

#+++++++++++
# Model 1 
#+++++++++++
immMNL_imm1 <- -2 * (immMNL$logLik + imm1$Loglik)
p.val_immMNL_imm1 <- pchisq(immMNL_imm1, df = 1, lower.tail = FALSE)

taxMNL_tax1 <- -2 * (taxMNL$logLik + tax1$Loglik)
p.val_taxMNL_tax1 <- pchisq(taxMNL_tax1, df = 1, lower.tail = FALSE)

ccMNL_cc1 <- -2 * (ccMNL$logLik + cc1$Loglik)
p.val_ccMNL_cc1 <- pchisq(ccMNL_cc1, df = 1, lower.tail = FALSE)


#+++++++++++
# Model 2 
#+++++++++++
taxMNL2_tax2 <- -2 * (taxMNL2$logLik + tax2$Loglik)
p.val_taxMNL2_tax2 <- pchisq(taxMNL2_tax2, df = 3, lower.tail = FALSE)

#----------------------------------------------------------
# Table 4: Likelihood Ratio Tests: Standard MNL vs. GHMNL
#----------------------------------------------------------
lrtest_stat  <-round(rbind(immMNL_imm1, taxMNL_tax1, ccMNL_cc1,taxMNL2_tax2), digits=3)
lrtest_p.val <-round(rbind(p.val_immMNL_imm1, p.val_taxMNL_tax1, p.val_ccMNL_cc1, p.val_taxMNL2_tax2), digits=3)


Table4 <- cbind(lrtest_stat, lrtest_p.val)
colnames(Table4)<- c("Test Statistic", "p-value")
rownames(Table4)<- c("Table 1: Immigration Issue: MNL vs. GHMNL", 
                            "Table 1: Tax Issue: MNL vs. GHMNL", 
                            "Table 1: Climate Change Issue: MNL vs. GHMNL", 
                            "Table 2: Full Model: MNL vs. GHMNL")
Table4


#==================================================================
# S-MNL Models
#==================================================================
dat_imm2 <- mlogit.data(dat_imm, choice = "chosen", shape = "long", alt.levels =c("1", "2", "3", "4", "5", "6"))
dat_tax2 <- mlogit.data(dat_tax, choice = "chosen", shape = "long", alt.levels =c("1", "2", "3", "4", "5", "6"))
dat_cc2  <- mlogit.data(dat_cc,  choice = "chosen", shape = "long", alt.levels =c("1", "2", "3", "4", "5", "6"))

# Set seed for reproducibility
set.seed(294271)

##############################
# S-MNL (without covariates)
##############################

# Model 1:
immSMNL1a <- gmnl(chosen ~ imm | 1, data = dat_imm2, model = "smnl", R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(immSMNL1a)

taxSMNL1a <- gmnl(chosen ~ tax | 1, data = dat_tax2, model = "smnl", R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(taxSMNL1a)

ccSMNL1a <- gmnl(chosen ~ cc | 1, data = dat_cc2, model = "smnl", R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(ccSMNL1a)


# Model 2: 
taxSMNL2a <- gmnl(chosen ~ tax | kath + gender + union |0, data = dat_tax2, model = "smnl", R = 500, 
                  notscale = c(1, 1, 1, 1, 1,0, rep(0, 15) ))
                                                                                                                
summary(taxSMNL2a)


##############################
# S-MNL (with covariates)
##############################

# Model 1:
immSMNL1 <- gmnl(chosen ~ imm | 1| 0 |0 | pimm -1, data = dat_imm2, model = "smnl", R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(immSMNL1)

taxSMNL1 <- gmnl(chosen ~ tax | 1| 0 |0 | ptax -1, data = dat_tax2, model = "smnl", R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(taxSMNL1)

ccSMNL1 <- gmnl(chosen ~ cc | 1| 0 |0 | pcc -1, data = dat_cc2, model = "smnl", R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(ccSMNL1)

# Model 2: 
taxSMNL2 <- gmnl(chosen ~ tax | kath + gender + union| 0 |0 |edubin + gender + west -1, 
                  data = dat_tax2, model = "smnl", R = 500, notscale = c(1, 1, 1, 1, 1,0,  rep(0, 15) ))
                                                                                                                
summary(taxSMNL2)

#############################################################################################################


#==================================================================
# MXL Models
#==================================================================

# Set seed for reproducibility
set.seed(294271)

##############################
# MXL (without covariates)
##############################

# Model 1:
immMXL1a<- gmnl(chosen ~ imm | 1 | 0 , data = dat_imm2, model = "mixl", ranp = c(imm = "ln"), R = 500)
summary(immMXL1a)

taxMXL1a<- gmnl(chosen ~ tax | 1 | 0 , data = dat_tax2, model = "mixl", ranp = c(tax = "ln"), R = 500)
summary(taxMXL1a)

ccMXL1a<- gmnl(chosen ~ cc | 1 | 0 , data = dat_cc2, model = "mixl", ranp = c(cc = "ln"), R = 500)
summary(ccMXL1a)


# Model 2
taxMXL2a<- gmnl(chosen ~ tax | kath + gender + union | 0 , data = dat_tax2, model = "mixl", ranp = c(tax = "ln"),  R = 500)
summary(taxMXL2a)


##############################
# MXL (with covariates)
##############################

# Model 1:
immMXL1<- gmnl(chosen ~ imm | 1 | 0 | pimm  -1, data = dat_imm2, model = "mixl", ranp = c(imm = "ln"), mvar = list(imm = c("pimm")), R = 500)
summary(immMXL1)

taxMXL1<- gmnl(chosen ~ tax | 1 | 0 | ptax  -1, data = dat_tax2, model = "mixl", ranp = c(tax = "ln"), mvar = list(tax = c("ptax")), R = 500)
summary(taxMXL1)

ccMXL1<- gmnl(chosen ~ cc | 1 | 0 | pcc  -1, data = dat_cc2, model = "mixl", ranp = c(cc = "ln"), mvar = list(cc = c("pcc")), R = 500)
summary(ccMXL1)

# Model 2
taxMXL2<- gmnl(chosen ~ tax | kath + gender + union | 0 | edubin + gender + west -1, data = dat_tax2, model = "mixl", 
               ranp = c(tax = "ln"), mvar = list(tax = c("edubin","gender", "west")), R = 500)
summary(taxMXL2)


#==================================================================
# G-MNL Models
#==================================================================

# Set seed for reproducibility
set.seed(294271)

##############################
# G-MNL (without covariates)
##############################

# Model 1
immGMNL1a <- gmnl(chosen ~ imm | 1| 0 |0 , data = dat_imm2, model = "gmnl", ranp = c(imm = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(immGMNL1a)

taxGMNL1a <- gmnl(chosen ~ tax | 1| 0 |0 , data = dat_tax2, model = "gmnl", ranp = c(tax = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(taxGMNL1a)

ccGMNL1a <- gmnl(chosen ~ cc | 1| 0 |0 , data = dat_cc2, model = "gmnl", ranp = c(cc = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(ccGMNL1a)

# Model 2
taxGMNL2a <- gmnl(chosen ~ tax |  kath + gender + union| 0 |0, data = dat_tax2, model = "gmnl", ranp = c(tax = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0, rep(0, 15) ))
summary(taxGMNL2a)


##############################
# G-MNL (with covariates)
##############################

# Model 1
immGMNL1 <- gmnl(chosen ~ imm | 1| 0 |0 |pimm  -1, data = dat_imm2, model = "gmnl", ranp = c(imm = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(immGMNL1)

taxGMNL1 <- gmnl(chosen ~ tax | 1| 0 |0 |ptax  -1, data = dat_tax2, model = "gmnl", ranp = c(tax = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(taxGMNL1)

ccGMNL1 <- gmnl(chosen ~ cc | 1| 0 |0 |pcc  -1, data = dat_cc2, model = "gmnl", ranp = c(cc = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0))
summary(ccGMNL1)


# Model 2
taxGMNL2 <- gmnl(chosen ~ tax |  kath + gender + union| 0 |0 |edubin + gender + west -1, data = dat_tax2, model = "gmnl", 
                 ranp = c(tax = "ln"), R = 500, notscale = c(1, 1, 1, 1, 1,0, rep(0, 15) ))
summary(taxGMNL2)


#------------------------------------------
# Table 5: Empirical Model Comparisons
#------------------------------------------

IMM1_compare <- matrix(NA, 8, ncol = 3)
colnames(IMM1_compare) <- c("LogL",  "AIC", "BIC")
rownames(IMM1_compare) <- c("MNL", "MXL (without covariates)", "MXL", "S-MNL (without covariates)", "S-MNL", "G-MNL (without covariates)", "G-MNL", "GHMNL")
IMM1_compare[1,] <- cbind(immMNL$logLik, AIC(immMNL), AIC(immMNL, k = log(1807)))
IMM1_compare[2,] <- cbind(immSMNL1a$logLik$maximum, AIC(immMXL1a), BIC(immMXL1a))
IMM1_compare[3,] <- cbind(immMXL1$logLik$maximum, AIC(immMXL1), BIC(immMXL1))
IMM1_compare[4,] <- cbind(immSMNL1a$logLik$maximum, AIC(immSMNL1a), BIC(immSMNL1a))
IMM1_compare[5,] <- cbind(immSMNL1$logLik$maximum, AIC(immSMNL1), BIC(immSMNL1))
IMM1_compare[6,] <- cbind(immGMNL1a$logLik$maximum, AIC(immGMNL1a), BIC(immGMNL1a))
IMM1_compare[7,] <- cbind(immGMNL1$logLik$maximum, AIC(immGMNL1), BIC(immGMNL1))
IMM1_compare[8,] <- cbind(-imm1$Loglik, (-2* (-imm1$Loglik) + 2*7), (-2* (-imm1$Loglik) + log(1807)*7))
IMM1_compare


TAX1_compare <- matrix(NA, 8, ncol = 3)
colnames(TAX1_compare) <- c("LogL", "AIC", "BIC")
rownames(TAX1_compare) <- c("MNL", "MXL (without covariates)", "MXL", "S-MNL (without covariates)", "S-MNL", "G-MNL (without covariates)", "G-MNL", "GHMNL")
TAX1_compare[1,] <- cbind(taxMNL$logLik, AIC(taxMNL), AIC(taxMNL, k = log(1251)))
TAX1_compare[2,] <- cbind(taxMXL1a$logLik$maximum, AIC(taxMXL1a), BIC(taxMXL1a))
TAX1_compare[3,] <- cbind(taxMXL1$logLik$maximum, AIC(taxMXL1), BIC(taxMXL1))
TAX1_compare[4,] <- cbind(taxSMNL1a$logLik$maximum, AIC(taxSMNL1a), BIC(taxSMNL1a))
TAX1_compare[5,] <- cbind(taxSMNL1$logLik$maximum, AIC(taxSMNL1), BIC(taxSMNL1))
TAX1_compare[6,] <- cbind(taxGMNL1a$logLik$maximum, AIC(taxGMNL1a), BIC(taxGMNL1a))
TAX1_compare[7,] <- cbind(taxGMNL1$logLik$maximum, AIC(taxGMNL1), BIC(taxGMNL1))
TAX1_compare[8,] <- cbind(-tax1$Loglik, (-2* (-tax1$Loglik) + 2*7), (-2* (-tax1$Loglik) + log(1251)*7))
TAX1_compare


CC1_compare <- matrix(NA, 8, ncol = 3)
colnames(CC1_compare) <- c("LogL", "AIC", "BIC")
rownames(CC1_compare) <- c("MNL", "MXL (without covariates)", "MXL", "S-MNL (without covariates)", "S-MNL", "G-MNL (without covariates)", "G-MNL", "GHMNL")
CC1_compare[1,] <- cbind(ccMNL$logLik, AIC(ccMNL), AIC(ccMNL, k = log(1807)))
CC1_compare[2,] <- cbind(ccMXL1a$logLik$maximum, AIC(ccMXL1a), BIC(ccMXL1a))
CC1_compare[3,] <- cbind(ccMXL1$logLik$maximum, AIC(ccMXL1), BIC(ccMXL1))
CC1_compare[4,] <- cbind(ccSMNL1a$logLik$maximum, AIC(ccSMNL1a), BIC(ccSMNL1a))
CC1_compare[5,] <- cbind(ccSMNL1$logLik$maximum, AIC(ccSMNL1), BIC(ccSMNL1))
CC1_compare[6,] <- cbind(ccGMNL1a$logLik$maximum, AIC(ccGMNL1a), BIC(ccGMNL1a))
CC1_compare[7,] <- cbind(ccGMNL1$logLik$maximum, AIC(ccGMNL1), BIC(ccGMNL1))
CC1_compare[8,] <- cbind(-cc1$Loglik, (-2* (-cc1$Loglik) + 2*7), (-2* (-cc1$Loglik) + log(1193)*7))
CC1_compare

TAX2_compare <- matrix(NA, 8, ncol = 3)
colnames(TAX2_compare) <- c("LogL", "AIC", "BIC")
rownames(TAX2_compare) <- c("MNL", "MXL (without covariates)", "MXL", "S-MNL (without covariates)", "S-MNL", "G-MNL (without covariates)", "G-MNL", "GHMNL")
TAX2_compare[1,] <- cbind(taxMNL2$logLik, AIC(taxMNL2), AIC(taxMNL2, k = log(1251)))
TAX2_compare[2,] <- cbind(taxMXL2a$logLik$maximum, AIC(taxMXL2a), BIC(taxMXL2a))
TAX2_compare[3,] <- cbind(taxMXL2$logLik$maximum, AIC(taxMXL2), BIC(taxMXL2))
TAX2_compare[4,] <- cbind(taxSMNL2a$logLik$maximum, AIC(taxSMNL2a), BIC(taxSMNL2a))
TAX2_compare[5,] <- cbind(taxSMNL2$logLik$maximum, AIC(taxSMNL2), BIC(taxSMNL2))
TAX2_compare[6,] <- cbind(taxGMNL2a$logLik$maximum, AIC(taxGMNL2a), BIC(taxGMNL2a))
TAX2_compare[7,] <- cbind(taxGMNL2$logLik$maximum, AIC(taxGMNL2), BIC(taxGMNL2))
TAX2_compare[8,] <- cbind(-tax2$Loglik, (-2* (-tax2$Loglik) + 2*7), (-2* (-tax2$Loglik) + log(1251)*7))
TAX2_compare

FullModel<-matrix(NA, 1,3)
rownames(FullModel)<-c("Full Model: Tax Issue")
FullModel


Table5 <-round(rbind(IssueImm,IMM1_compare, IssueTax,TAX1_compare, IssueCC, CC1_compare, FullModel, TAX2_compare), digits=3)
Table5
